Finite differences schemes

To differentiate some field(s) in one direction using finite differences we provide a generic stencil generator facility StencilGenerator.

If tab is a numpy array representing the discrete values of a scalar field \(\rho\) on a grid, then to compute

\[\begin{eqnarray} result &=& \frac{\partial \rho}{\partial y} \end{eqnarray}\]

the basic usage of such schemes is a two step process:

  1. Generate a stencil given dimension, derivative and order.

Example with 1st derivative order 2 in 1D:

>>> from hysop.numerics.stencil.stencil_generator import StencilGenerator
>>> sg = StencilGenerator()
>>> st = sg.configure(dim=1, derivative=1, order=2).generate_exact_stencil(origin=1)
>>> print(f"Order {st.order}, {st.shape[0]} points scheme with coefficients :{st.coeffs}")
Order 3, 3 points scheme with coefficients :[-1/2 0 1/2]

Example with 2nd derivative order 2 in 2D:

>>> sg = StencilGenerator()
>>> st = sg.configure(dim=2, derivative=2, order=2).generate_exact_stencil(origin=1)
>>> print(f"Order {st.order} {st.shape} points scheme with coefficients :{st.coeffs}")
Order [3 3] (3, 3) points scheme with coefficients :[[0 dx**(-2) 0]
 [dy**(-2) (-2*dx**2 - 2*dy**2)/(dx**2*dy**2) dy**(-2)]
 [0 dx**(-2) 0]]

Stencil has been generated with dx and dy symbols. Real values can be replaced instead:

>>> dx,dy=st.dx
>>> st.replace_symbols({dx: 1, dy: 1})
>>> print(st.coeffs)
[[0 1 0]
 [1 -4 1]
 [0 1 0]]

  1. Apply stencil to array using apply method

>>> import numpy as np
>>> sg = StencilGenerator()
>>> st = sg.configure(dim=1, derivative=1, order=2).generate_exact_stencil(origin=1)
>>> st.replace_symbols({st.dx: 0.1})
>>> st.apply(np.linspace(0,1,11))[1:-1]
array([1., 1., 1., 1., 1., 1., 1., 1., 1.])

A few important remarks

  • Mind that apply first argument is a numpy array not a field. Usually, you will have to use some_field.data[i] as tab, some_field being a DiscreteField and i the component number.

  • step is the space discretization step size in each direction, i.e. a list or numpy array with d values, d being the dimension of the domain. A common case is code:

    step = topo.mesh.space_step

  • indices represent the set of points on which the scheme must be applied. This is usually a list of slices, for example, code:

    indices = topo.mesh.compute_index

  • result must be a predefined numpy array of the ‘right’ size, here the same size/shape as tab.

  • To optimize memory usage and computational time, it’s possible to reduce the size of the output and/or to apply the scheme on a subset of the domain. scheme.indices represents the indices (in topo.mesh) of points for which the derivative will be computed. scheme.output_indices are positions in result where the resulting derivative will be saved. All available possibilities are summarized through the examples below.

  • The size of the ghost layer depends on the scheme but is not checked! You must ensure that topo.ghosts() >= scheme.ghosts_layer_size.

  • For some schemes a work array may be provided during call. If so, it must be a numpy array of the same size as result. It’s shape is not really important, since during call work will be reshaped to be of the same shape as result. This allows us to provide 1D array that may be shared between different operators/methods whatever the internal required shapes are.

  • Notice that most of the time, finite-differences are defined as internal methods of operators and work arrays management, indices list or ghosts layers are set/checked internally.

Default case : apply scheme on all the domain with a full-size result

This is the standard case where result and input tab have the same shape. Computation is performed on all grid points, excluding ghosts.

>>> import numpy as np
>>> from hysop import Box, Field
>>> from hysop.topology.cartesian_topology import CartesianTopology
>>> from hysop.tools.parameters import CartesianDiscretization
>>> # --- Domain, discretization and topology ---
>>> box = Box(length=[1., 1.])
>>> d2d = CartesianDiscretization([33, 33], [1, 1], default_boundaries=True)
>>> topo = CartesianTopology(box, d2d)
>>> # --- Scalar field on the domain ---
>>> rho = Field(box, name='rho')
>>> # discretize and initialize rho
>>> rd = rho.discretize(topo)
>>> _ = rd.randomize()
>>> # --- Grid parameters required for FD scheme ---
>>> # space step
>>> step = topo.mesh.space_step
>>> # field resolution on the grid defined by topo
>>> shape = topo.mesh.local_resolution
>>> # create a new array to save results
>>> result = np.zeros(shape)
>>> # --- Create FD scheme ---
>>> sg = StencilGenerator()
>>> st = sg.configure(dim=2, derivative=1, order=2).generate_exact_stencil(origin=1)
>>> st.replace_symbols({st.dx[0]: step[0], st.dx[1]: step[1]})
>>> # --- Apply FD scheme ---
>>> # compute first derivative of scal on all points (except ghosts)
>>> # of the 'topo' discretization of the domain
>>> result = st(rd.sbuffer,result, (0,1))

Notice that ghost points of result are not updated!

Apply scheme on all the domain with a reduced result

Suppose that you want to compute a derivative on all the domain but use non-standard shape array to save the result. This may be useful if you do not want to allocate memory for ghost points in result.

>>> shape = topo.mesh.compute_resolution
>>> result = np.zeros(shape)
>>> # Get 'computational' points, i.e. grid points excluding ghosts
>>> ic = topo.mesh.local_compute_slices
>>> result = st(rd.sbuffer,result, (0,1), iview=ic)